Numerical Implementation of Generalized Unitarity 



Petros Draggioti^ll 

University of Athens, Physics Department, Nuclear & Particle Physics Section, Athens, Greece 



Abstract 

Wc present a numerical implementation of generalized unitarity. Wc will show that we 
are able to obtain the box coefficients of any 1 loop gluonic amplitudes for an arbitrary 
helicity configuration and for any number of external particles. 

1 Introduction 

One loop amplitudes is an important and very complicated part of a scattering cross section 
calculation in high energy physics. Even for a small number of external particles the number of 
Feynman diagrams grows really fast, when one adds just one more leg. During the last years, 
there has been a lot of progress in understanding the complex structure of 1 loop amplitudes 
and new methods have helped us to work around Feynman diagrams in getting to the result. 
These methods depend on the complex plane structure of an amplitude when we continue all 
momenta to become complex. 

The particular method that we focus, generalized unitarity, is a clever method of obtaining 
the coefficients of the master integrals, that constitute the diverging part of an amplitude. 
Based on cutting the propagators in the loop and obtaining a particular coefficient as product of 
tree amplitudes, the method lends itself naturally for an algorithmic implementation. We have 
developed just such an implementation, which (as a first step) computes the box coefficients of 
any 1 loop amplitude for any helicity configuration of the external particle. Since, implicitely, 
the method relies on the calculation of tree amplitudes with various helicity configurations, we 
have also numerically implemented the BCFW recursion equation. 

2 The BCFW recursion equations 

In this section we briefiy present the BCFW recursion and its numerical implementation. But 
first let us agree on some definitions and convention that we use throughout this work. 

2.1 Spinor definitions and conventions 

It is a well established fact how a massless four-vector can be written as matrix using the Pauli 
matrices a^: 

Pi^ Pab = P^^'^AB 
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Being a massless vector, the matrix p^^ has a determinant zero det{p) = 0. It is known from 
linear algebra that a matrix with zero determinant can be written as a product of two columns 
which we call spinors: 

Pab = ^a\b (2) 
We write these in a convienient bra-ket form as follows: 

^A{Pi) V) , A^(Pj) {i\ (3) 
For definite helicities we can write these spinors as solutions of the Weyl equation: 

\i±) = u±{pi) {i±\ = u±{pi) (4) 

With these definitions we can define inner products with spinors. There are two versions the 
'angle' product and the 'bracket'. For real momenta these are conjugate of each other. For 
complex momenta this may no longer hold. 

(ij) = =u^{pi)u+{pj) [ij] = + = u+{pi)u^{pj) (5) 

From the definitions of the spinors we can see that products like {i + or {i — vanish. 
We can use these definitions to construct more complicated objects, with momenta sandwiched 
between spinors: 

(i + \k\j+) = {i + = [ik\{k3) {i - \k\j-) = {i- W\j-) = {ik)[kj] (6) 

with a generalization for n momenta, 

{i + \k1k2 ■ ■ ■ kn\j+) = {i + Wi--- ¥n\j+) = [iki]{kik2)[k2k3] ■ ■ ■ {kn-ikn) (7) 
We will also use the identities: 

n 

{i + \ki + k2 + --- + kn, + --- + kn \j+) = [im] {mj) (8) 

m=l 

n 

{i-\ki + k2 + --- + km + --- + kn\j-) = {im)[mj] (9) 

m=l 

Finally, as you may have already noticed, we will use the index of the momentum as the symbol 
for the spinor, when there is no confusion. 

2.2 Analytic part 

For a LO amplitude, and concentrating on gluons only for simplicity, the color information can 
be factorized from the kinematical part as follows: 

A*r%{ki, Xi}) = g""-^ rr(rMi)TM2) ...TM"))A^^e^(a(l^i),...,(T(n^")) (10) 
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where are the SU{N) matrices in the adjoint respresentation. The sum is over ah non-cyclic 
permutations of the external legs. The part of the factorization A!^'^^{a{l^^), . . . ,c7(n'*'")) now 
depends only on the external momenta, and the polarization vectors and is in a color ordered 
form. This has the nice property that all poles that appear in it are made of adjacent momenta, 
■si2,'S23 etc. An almost identical result holds for quarks as well. 

Some years ago , Britto, Cachazo and Feng [1], showed that the kinematical part of the 
amplitude, the so called partial amplitude, has a nice factorization property that satisfies the 
recursion equation ( thereafter known as the BCFW recursion): 

71-2 

<-(l,2, ...,n)= J]47^(i,2, . . . , -K-^^)-^A^lU,{Kl„ . . . ,n - l,n) (11) 

h=±k=2 l.fc 

The color ordered amplitude is split into two terms, consisting of lower point ordered amplitudes, 
joint by a propagator Ki^k = Pi + P2 + ■ ■ ■ + Vk- With an appropriate complex shift of the 
momenta the propagator can be put on shell. The simplest way to do that is by shifting the first 
and the last leg (denoted by a hat in the recursion), by a fixed amount with a complex constant 
z: 

1^1 = 1 + z^n , 1^1 

n ^ n , n ^ n = h — Zkl (12) 

The z is chosen in such a way so that the shifted propagator is massless: 

= ~ 7 17^' n \ (13) 

(n - ||ri,A;|l-) 

Klk = + ^kulf = Kl, + Zk{n- ||ri,fc|l-) = (14) 

The sums are over all possible distributions of the external momenta into the two groups, 
keeping the special, shifted momenta in distinct groups, and over the helicities of the common 
propagator. The recursion was proven by Britto et al [2] using the analytic properties of the 
amplitude and by Lazopoulos et al [3] using standard Feynman diagrams analysis. 

The advantage of the BCFW recursion is that only a few terms are needed to compute a 
partial amplitude, compared to hundrends of Feynman diagrams using standard field theory 
techniques. This signals major cancellations between graphs that are automatically taken care 
of using the recursion. Thus, we arrive at relatively compact expressions for the amplitudes. 



2.3 Numerical implementation 

We have implemented the BCFW recursion, in a numerical F0RTRAN90 code, that computes 
the partial amplitudes for QCD, including both quarks and gluons, for an arbitrary number 
of external legs. The only inputs in the code is the number of legs n, the QCD process, and 
the helicity configuration. We have compared the results, both for MHV and NMHV helicity 
configurations with known results from the literature and with HELAC 1^. We have also 
produced some new numerical results. Some results are shown in the following tables for 8 and 
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9 point amplitudes with and without quarks. We have used the following momenta values for 
the 8-point amplitude: 



Pi = 


(50.0,0.0,0.0,50.0) 




P2 = 


(50.0,0.0,0.0,-50.0) 




P3 = 


(-11.1703767834197,- 


8.03205150805747, 6.45958520425105, -4.30548778241159) 


P4 = 


(-29.3374489041782, - 


0.341091815726772, 5.38045260205076, 28.8378274905652) 


P5 = 


(-8.93691541288138, - 


7.79003576692047, 2.54732281477208, -3.56299681791019) 


P6 = 


(-4.49061013802485, - 


2.19138756179328, 0.744707846345424, 3.84822169181309) 


P7 = 


(-9.04125786194674,4 


.56116168271723, 7.76548604819915, 0.798357229284189) 


P8 = 


(-37.0233908995491, 13.7934049697808, -22.8975545156185, -25.6159218113407) 



Helicity Configuration 


Partial Amplitude 


(l-,2-,3-,4-,5+,6+,7+,8+) 


1.143995504515313 10"* + 2.452584716117081 IQ-^i 


(l-,2+,3+,4-,5-,6+,7+,8~) 


-3.892700175442802 10"** - 7.138278322124297 IQ-^i 


(l-,2~,3-,4+,5+,6+,7+,8+) 


2.973212107227238 lO"® + 1.495525686984027 lO^^i 


(l-,2-,3+,4+,5-,6-,7+,8+) 


-8.620198212038414 10"** + 8.089172271501075 IQ-^i 


(l+,2+,3+,4-,5+,6-,7+,8-) 


1.899386258601871 lO^^ - 9.145607147903538 IQ-^i 


(lq,2+,3+,4+,5+,6-,7-,8-) 


-1.812332579250756 10"^ + 1.534899358766671 IQ-'^i 


(lq,2±,3-,4+,5-,6+,7-,8+) 


-5.724656620620708 lO^^ - 8.875798503634412 10-^°i 



Table 1. Results for 8-point partial amplitudes. Quarks and antiquarks are denoted as subscripts q and q in 

the relevant legs 



The following momenta configuration was used for the 9-point amplitudes: 



Pi = 


(50.0,0.0,0.0,50.0) 


P2 = 


(50.0,0.0,0.0,-50.0) 


P3 = 


(-8.62171654444322, -6.90534943082250, 4.59215344257906, -2.35844698282037) 


PA = 


(-27.8413520112470, -0.643296108254438, 0.840299169363624, 27.8212319863308) 


P5 = 


(-7.13369828021707, -6.67842455614503, 1.44431114936003, -2.04994192719828) 


P6 = 


(-4.22769845708132, -1.90930333047348, 7.776353148874293 10"^ 3.77119979176668) 


P7 = 


(-7.10934176699029, 3.77076042627538, 5.74622866522395, 1.81795552598701) 


P8 = 


(-30.9459371413082, 11.4021700192282, -22.0566594617972, -18.4701195933519) 


P9 = 


(-14.1202557987129, 0.963442980191897, 9.35590350378182, -10.5318788007139) 
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Helicity Configuration 


Partial Amplitude 


(l-,2-,3-,4-,5+,6+,7+,8+,9+) 


-1.195041826437500 10"*^ - 4.267905334966595 lO^^^i 


(l-,2+,3+,4-,5-,6+,7+,8-,9-) 


9.960596678335883 10"^ + 2.632108576287432 lO^^^i 


(l-,2-,3-,4+,5+,6+,7+,8+,9+) 


-2.628610959245734 lO^^" - 3.014249899217379 IQ-^i 


(l-,2-,3+,4+,5-,6-,7+,8+,9-) 


-8.470520978313771 lO"'^ - 1.528306470088502 lO^^i 


(l+,2+,3+,4-,5+,6-,7+,8-,9+) 


-1.935502160225642 10"^ + 7.283711249910175 10~^i 


(lq,2t,3+,4+,5+,6-,7-,8-,9-) 


-4.751617821378328 10"^ - 5.102361545628036 10"^^ 


(lq,2^,3-,4+,5-,6+,7-,8+,9-) 


-4.359890532350093 10"^ + 1.991698153552900 lO^^i 



Table 2. Results for 9-point partial amplitudes. Notations are the same as before 



3 The next step: NLO 

The color factorization structure for the LO amplitudes, persists for NLO amplitudes only 
slightly more involved: 

n 

n/2+l 

+ 5" 5] E Gn;c{<^)An.,c{a{l^'),...,a{n^")) (15) 

C=2 a€:Sn/Sn;c 

where Gn;i{cr) is the leading color structure: 

Gn-i{(T) = iVerr(T''-(i)T"-(2)---T"-(")) (16) 

and Gn-c{(^) is the subleading color structure: 

Gn;c{(^) = Tr {T^-w . . . T'^<.(e-i) ) Tr {T"-^-) ■ ■ ■ T"-(") ) (17) 

The colored ordered amplitudes of the subleading terms, ^ri;c(c), can be written as linear com- 
binations of permumations of the leading color amplitudes An-i[a)^ so the latter are called 
primitive amplitudes. For the rest of this paper we focus on these primitive amplitudes. 

By restricting ourselves to cyclicaly ordered primitive amplitudes we reduce the labour of 
computing hundrends or thousands of Feynman diagrams to a reduced set of diagrams. Now 
each of those diagrams is an integral over the loop momentum. The integrands of those integrals 
are made of tensors or vectors of the external momenta and the loop momenta. It would an 
enormous simplification if we could further reduce these to a basic set of integrals, so that every 
diagram could be written as a combination in this basis. It turns out that such a reduction 
is possible, either using standard Passarino-Veltman techniques [5], or the more recent Ossola- 
Papadopoulos-Pittau [6] reduction method. Thus all loop integrals can be brought down to a 
basis consisting of up to 4 propagators: 

T fT Tim j2m jSm jlm j2m e j2m h j3m r4mi 

-L — U2,^3 '-'3 '-'S '-'4 '-'4 )-'4 ;-'4 ; -'4 I UoJ 
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As a result any 1 loop amplitude can be written as a linear combination in this basis, with 
algebraic coefficients. 



(19) 

j 

where TZn are the rational terms. More explicitly we can write: 

^i-ioop ^ ^ ^.^(^2) ^ J2 c,,CiKf, K]) + d^JkD{Kf, K],Kl) + 7^„ (20) 

i ij ijk 

where B, C and D are the bubble, triangle and box integrals respectively. Analytic expressions 
and singularity structure for these integrals can be found for example in [7] 

4 Generalized Unitarity 

Unitarity of the scattering matrix in field theory implies the conservation of probability. On 
the amplitude level unitarity says that cutting a loop gives the discontinuity in the scattering 
amplitude. By cutting we mean putting a propagator on shell, which amounts in the replacement: 

27r5W(/) (21) 



+ ie 



At the diagrammatic level this the well known Cutkosky rule. In general, at the amplitude 
level, cutting an amplitude in a given channel isolates those integrals that have a discontinuity 
across that channel. In other words cutting gives a linear combination of master integrals with 
algebraic coefficients. If we wish to isolate one single integral we must go a bit further and insist 
more propagators go on shell. This goes under the name of generalized unitarity |S] . For the box 
integrals in particular this is quite easy to do. Cutting four propagators, isolates a single box 
integral. The cut breaks the loop integral in four tree level amplitudes (one at each corner of 
the box) and the coefficient of the box integral is simply the product of those tree amplitudes: 



• \ 


^ 1 


/> 




1 
1 


1 I 


1 













Fig.l Cutting 4 propagators, gives the Box coefficient, as a product of the tree level 
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a=l h=± 

x4-^(4^ {fc} , -C-)A'r{&-,{m} , -l-^^) (22) 

where {?} , {i} etc., denotes the collection of momenta in that vertex of the box. One also 
averages over the two solutions of the quadratic equation one obtains by putting four propagators 
on shell. Finally there is a sum over the helicities of the propagators. 

Simple example ^(l"'', 2+, 3~, 4~) We have to satisfy the on-shell conditions: l\ = 0,/! = 
0,^3 = 0,^4 = 0. Choosing I2 to be : ^2 ~ '^(^ ~ It'^I^— ) ) where ^ is a normalization, satisfies 
the first three. The last will determine ^. The product of the trees in the corners of the box is: 



1 [1/2]^ (^44)3 (/3/2)=^ [hh] 



(23) 



2[hh][hl\ {4h){hh) {h2){2h) [hmU] 
We can combine the numerator into one factor: 

[l/2](?2/3) 1/3/4] (/44) = (1 + 1/2/3/414+) (24) 

Using the identities Ijil/i = + ^ = 10 and I/^IM = = /^(/i + '^^ = l/^W this can be 
simplified: 

(1 + I/2/3/4I4+) = (1 + |/223|4+) = (1 + |/2|2+)[23](34) = si2[23](34) (25) 

Combining appropriately factors in the denominator, we can simplify the expression: 

[3/4](/4/i)[/i/2](/22) = (3 + I/4/1/2I2+) = [34](41)(1 + I/2I2+) = [34](41)si2 (26) 

(4/i)[/il] = (4-|/i|l-) = (41)[21] (27) 
(2/3)[/33] = (2- I/3I3-) = (21)[23] (28) 
Putting together numerator and denominator and using conservation of momenta we finally get: 

4:2:3:4 = ^Si2S234'''(l"' , 2+ , 3", 4") (29) 

5 Numerical Implementation 

We have implemented generalized unitarity in a numerical FORTRAN 95 code. For the time being 
it has been implemented for box coefficients only and for purely gluonic amplitudes (no quarks 
in the loop). The evaluation of the box coefficient goes through the following steps: 

• Lists all possible groupings for the external momenta: The external momenta are 
grouped in 4 groups (the corners of the box), in all possible, distinct ways. Sums of these 
momenta make up the channel that the cut is computed. 
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• Solutions of the onshell loop momenta condition: The system of equations: if = 
0, £2 = 0, £| = 0, £| = for the cut loop momenta is solved and the two solutions, that we 
have to average over are evaluated 

• The tree amplitudes are computed: For each corner of the box we compute the 
tree amplitude defined by the external momenta in that corner, the solution for the loop 
momentum that we just obtained and the particular helicity configuration of the grouping 
of the external momenta. The computation of the tree amplitudes is done automatically 
using the numerical code for the BCFW recursion that we discussed in a previous section. 

• Finally the coefficient for that particular cut is computed, using Eq. (|22|) . 

We have produced some results for two 1-loop processes, namely gg — > 8g and gg 9g 
for various choices of helicity configurations and for various channels (particular coefficients of 
a master integral). We have used the same momenta configuration as in the BCFW recursion 
section. We have compared our results with the ones in the literature and we have also produced 
some new numerical results in the case of 1-loop 9 gluon amplitudes. The results are shown in 
the following tables. 





- H h - + - + 


+ + ++ 


(1:2:3: 45678) 
(1 : 23 : 4 : 5678) 
(1:2: 345 : 678) 
(1 : 2 : 34 : 5678) 
(12 : 34 : 56 : 78) 


-0.452833202541873 + 1.83433531772824j 
-1.51120109392466410-5 _ 2.69252629629608410-^1 
-0.596872851993693 - 1.13690332085509* 
0.343744213463153 - 0.116177716850307i 
-6.08936506243673310-^ - 0.185567043177724* 


564.146943978796 + 341.450920347787* 
85.7803228327941 - 4.79606435167831* 
366.912336610440 + 98.0271230853793* 

0.0 + 0.0* 

0.0 + 0.0* 





-+-+-+++ 


- H h + H h 


(1:2:3: 45678) 
(1 : 23 : 4 : 5678) 
(1:2: 345 : 678) 
(1 : 2 : 34 : 5678) 
(12 : 34 : 56 : 78) 


-1.810002621799826 10-^ - 7.007553695597921 lO-^i 
-9.061239878511741 10-^ -6.182682488354774 IQ-'^i 
-6.198557976336275 10-^ - 3.145711006839414 lO"'^* 
2.057234610074154 10"^ + 5.080380495579885 10-^* 
0.0 + 0.0* 


1.282527393614916 IQ-^ - 1.544463085586650 10"^* 
-3.593566469301429 lO-'^ - 2.807114551959228 10-^* 
1.201810688581596 IQ-^ + 1.747770080527506 10-** 

0.0 + 0.0* 

0.0 + 0.0* 



Table 3. Results for 1-loop, 8-point primitive amplitudes. 



The code provides all the coefficients for a particular process at once. Typical times for the 
two processes shown in the tables are of the order of 1/lOth of a second on a 64-bit Intel Core 
Duo at 2GHz. 
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-H \ h- + - 


h + + + + 


(1:2:3: 456789) 
(1 : 23 : 4 : 56789) 
(1:2: 345 : 6789) 
(1 : 2 : 34 : 56789) 
(12 : 34 : 56 : 789) 


-0.291597431159727 + 0.263428066914340i 
6.917959935625083 10~^ - 2.862781163146907 lO^'^i 
1.83770214755652 + 1.25751379616759i 

-0.801669928958509 - 0.704349268560792i 
4.147713723224443 10^^ + 3.388946250063977 lO^^j 


-9.42379749599236 + 8.63580103044842i 
-4.424679500634447 10"^ - 1.06642829663490* 
-3.16966894074953 - 10.5119953247915i 
0.0 + O.Oi 
0.0 + O.Oi 





-H \ h + + - 


-H h + H h- 


(1:2:3: 456789) 
(1 : 23 : 4 : 56789) 
(1:2: 345 : 6789) 
(1 : 2 : 34 : 56789) 
(12 : 34 : 56 : 789) 


-2.382052305295742 10"^ - 1.320540695786225 lO^^i 
2.330322361904153 10"** - 2.184457500759925 lO"** 
1.731650767307649 lO'^ _ 1.381732810170124 lO'^j 
-9.470193967875743^0-=* + 8.652362025481619 lO'^i 
-8.488224296845019 lO"'"^ - 1.761227642104737 lO"** 


7.517076253761092 10"^ - 2.651249643255122 lO^^^i 
7.985298353983870 10"* - 1.068015483724008 lO-'^i 
1.377043202412121 lO"-* + 1.023698942546792 lO"'** 
-5.650383858159608 10""' - 1.479340069213934 ID'H 

3.287516684614978 10^** + 1.134952815474561 lO^'^i 



Table 4. Results for 1-loop, 9-point primitive amplitudes 



Triangle coefficients can be obtained in much the same way. Cutting 3 propagators does not 
immediately isolate a single integral but it gives a triangle integral plus a sum of box integrals [9] . 
This is because the particular box integrals share the same cuts with the triangles. The boxes 
have to be subtracted in a suitable way. This can naturally be implemented in our algorithm 
since box coefficients are computed already. For bubble integrals, cutting 2 propagators and 
subtracting the triangle and box contributions will isolate a single coefficient. Work is in progress 
to implement just this procedure. 

6 Summary and Outlook 

We presented numerical implementations for the BCFW recursion equations and 1-loop coeffi- 
cients of box integrals using generalized unitarity. In the case of the BCFW recursion equations 
we were able to provide, in a fast algorithmic way, multi-particle amplitudes both old and new, 
namely in the case of 9 parton Leading Order amplitude. For the implementation of generalized 
unitarity we are able to compute all box coefficients at once for an arbitrary number of particles 
and helicity configuration. Some new coefficients for the 1 loop 9 gluon amplitude were pre- 
sented. Work is in progress to include quarks in the picture (easy step since this is already done 
in the BCFW recursion code) and also complete the implementation by computing triangle and 
bubble coefficients. 
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